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Abstract 



For multiscale gas flows, kinetic-continuum hybrid method is usually used to balance the computa- 
tional accuracy and efficiency. However, the kinetic-continuum coupling is not straightforward since the 
coupled methods are based on different theoretical frameworks. In particular, it is not easy to recover 
the non-equilibrium information required by the kinetic method which is lost by the continuum model at 
the coupling interface. Therefore, we present a multiscale lattice Boltzmann (LB) method which deploys 
high-order LB models in highly rarefied flow regions and low-order ones in less rarefied regions. Since this 
multiscale approach is based on the same theoretical framework, the coupling precess becomes simple. The 



non-equilibrium information will not be lost at the interface as low-order LB models can also retain this in- 



q 

formation. The simulation results confirm that the present method can achieve model accuracy with reduced 
iy>j , computational cost. 

Si. 

i — i! 1 Introduction 

£> ■ Many engineering problems involve multi-scale gas flows, e.g. gas flows in micro/nano-fluidic devices 1 . Since 
JjL^ ■ the flow regions can be highly rarefied (non-equilibrium) , the conventional continuum theory becomes inappro- 
priate. The rarefaction order of gas flows can be classified by the non-dimensional Knudsen number, Kn, defined 



as the ratio of the mean free path and the device characteristic length scale. The Navier Stokes equations with 
no- velocity-slip wall boundary condition are only appropriate in the hydrodynamic regime where Kn < 0.001. 
In the slip flow regime (0.001 < Kn < 0.1), or the transition flow regime (0.1 < Kn < 10), it is necessary to use 
kinetic methods, e.g. direct simulation Monte Carlo method (DSMC) method, to describe gas flows. Although 
the kinetic methods including DSMC are able to simulate flow in the continuum or near-continuum regimes, the 



computational cost is often very expensive especially for low speed flows. Therefore, kinetic-continuum hybrid 
methods are naturally employed to deal with mixed flow regimes in typical microfluidic devices operating with a 
range of Knudsen numbers in different parts, i.e., numerically efficient continuum approach will be employed for 
the continuum regimes and kinetic approach for the rarefied regimes(see Refs. 0G2H1EJIE1II1IE1II2 EH HU HH US] 
and references therein). The two models are coupled by exchanging information over the "hand-shaking" region 
or across an interface. 

However, the kinetic-continuum coupling is not straightforward since two types of methods are based on 
different theoretical frameworks. While information transferring from the kinetic model to the continuum model 
is usually a well-defined process, the reverse process is more problematic 12 . It is difficult to recover non- 
equilibrium information lost by the continuum solvers which is required by the kinetic method. Although the 

1 



kinetic model can provide necessary information for the continuum model, it can be computationally expensive 
[3] . The statistical noise associated with the particle methods may also affect the accuracy and stability of the 
hybrid solver [T2]. To resolve these problems, we introduce a multiscale lattice Boltzmann (LB) method to use 
LB models which are based on a same theoretical framework. 

The LB method has been proven to be able to simulate hydrodynamic flows with only minimal number of 
discrete velocities (e.g., nine discrete velocities for a two-dimensional problem) [TU [T3 HH H3 [TS] . For continuum 
problems, its applicability is ensured by the Chapman-Enskog expansion. Due to its kinetic nature, the LB 
model has distinct advantages over other continuum computational methods, including easy implementation of 
multi-physical mechanisms and the boundary conditions for fluid/ wall interactions. 

The LB method may also offer a flexible framework for rarefied flows, which has recently been demonstrated 
extensively [HI [2HI [2H I2H [231 HH 123 [2H1 HZ1 [23 [23 1301 HZ]- It was shown that the key to capturing the 
rarefaction effects is to choose appropriate discrete velocity sets. Generally, a high-order LB model with larger 
discrete velocity set describes non-equilibrium effects better [23 123 123 OH] ■ In particular, high-order models 
with modest discrete velocity sets can already accurately capture non-equilibrium effects in rarefied flows over 
a range of Knudsen numbers [23 [23 123 EI] • 

Since the LB method offers a solution for simulating gas flows ranging from continuum to rarefied, we can 
introduce a multiscale method to couple models based on the same LB framework. This can be accomplished by 
employing higher-order LB models for non-equilibrium flow regions and lower-order LB models for hydrodynamic 
flow regions. Since the coupled LB models only differ in the chosen discrete velocities without loss of kinetic 
information at the coupling interface, the information exchange process can be simplified. In particular, non- 
equilibrium information can be retained in lower-order LB models, which resolves an intrinsic obstacle associated 
with kinetic-continuum hybrid methods. Meanwhile, the LB model can still reserve some advantages of particle 
method while eliminating the statistical noise. 

2 Lattice Boltzmann method 

2.1 Lattice Boltzmann equation 

Historically, the LB method was developed from the lattice gas cellular automata. The purpose was to mimic 
the Navier- Stokes dynamics. However, it was revealed that its applicability should not be limited to the 
hydrodynamic level [23E3E3E3E3H3I3I1EI112I]- There are different theoretical frameworks for LB models 
e.g., entropic LB models 37, 38, 39j.40j 41] . Here, to demonstrate the coupling scheme, we will adopt commonly- 
used LB models based on the Hermite expansion detailed in Refs.|26, 32, 33, 34, 35, 36 . However, the proposed 
coupling approach can be equally applied to different LB models. 
The original BGK equation is given as: 

% + S-vf + g-vd = - l U-f e9 ), (i) 

at [A 

where / denotes the distribution function, £ the phase velocity, p the pressure, g the body force and /i the 
gas viscosity. Using the well-known Chapman-Enskog expansion, the collision frequency can be represented by 
the ratio of pressure and gas viscosity, which is convenient to obtain the Knudsen number definition consistent 
with that of hydrodynamic models. Without losing generality, One can define the following non-dimensional 



variables: 



„ r . U y/M^t 
r = —,u= ,t = — , 





. _ ]o9_ £ _ g f=— (2) 

where u is the macroscopic velocity, R the gas constant, T the gas temperature, Tq the reference temperature, 
r the spatial position and Iq the characteristic length of the flow system. The symbol hat, which denotes a 
dimcnsionless value, will hereinafter be omitted. The Knudsen number can be defined by using macroscopic 
properties as below: 

Kn= - . (3) 

ph 

Based on these non-dimensional variables, the non-dimensional form of the BGK equation becomes 

f + *.V/ + ff -V,/ = -^(/-/-), (4) 

where the Maxwell distribution in D-dimensional Cartesian coordinates can be written as 

'-ii-ur 



r ~ (2nT)D/2 exp 



2T 



(5) 



For solving Eq.Q, the velocity space can be firstly discretized by projecting the distribution function onto 
a functional space spanned by the orthogonal Hermite basis [261142] : 



N 1 
/(r, 1 1) « f N (r, 1 1) = u{£) £ -a'"' (r, t) X {n) (0, (6) 



n=Q 



where x 1S the nth order Hermite polynomial, and w(£) is the weight function, which are given by 

f— Tl™ 

x ( " ) (0 = L 7 jrV"c(0, (7) 

««) = 7^757i^/ 2 . (8) 



(2tt) d / 2 



The coefficients o^™' are 
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[fX in) dS*[f W X in) dS (9) 






The equilibrium distribution should also be expanded as 
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f'^^Ei^to, (io) 



n=0 
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where the coefficient a e q for the equilibrium distribution is 



a 



<<i 



(n) _ / ceq In) 



f eq X w d£. (11) 



w a and £ a , a = 1, • • • , d, are the weights and abscissae of a Gauss- Hermite quadrature of degree > 2N respec- 
tively. Therefore, the Maxwell distribution is approximated by up to N Hermite polynomials. The body force 
term F(r, £,t) = g ■ Vj/ can also be approximated as[26l [35] 



N 1 

F(r,*,t) = w J2 T—Tr^^X^. (12) 

n=l ^ '" 

It was shown that Eq.Q with the first-order Hermite expansion is sufficient to capture the rarefaction effects 
for isothermal and incompressible flows [31] . On the other hand, the second order expansion has been proven to 
be able to model various Navier-Stokes level problems [T7J [T5]. Therefore, the second order approximation of 
the equilibrium distribution and the body force will be used hereinafter, as given below: 

r « ««M1 + i ■ u + i[(£ • uf - u 2 + (T l)(e D)]}, (13) 

F(r, £, t) » oj(i)p{g ■ * + (ff •€)(«•«)- ff • «}, (14) 

where T should be unity for isothermal problems and p is constant for incompressible problems. 

The discrete velocity set is revealed to be of upmost importance in determining model accuracy for rarefaction 
effects [3T]. For the Navier-Stokes level problems, several sets have been found to be applicable, e.g., the well 
known D2Q9j model for two-dimensional flows. To capture higher-order rarefaction effects, more discrete 
velocities are required. Some modest discrete velocity sets were shown (e.g., D2Q16 and D2Q36) to be able to 
capture non-equilibrium effects for flows over a broad range of Knudsen numbers [2"8l |3~T1 |2"TI |2"5] . Nevertheless, 
highly accurate discrete velocity set is required for the flows with large Knudsen number. Therefore, coupling 
high-order and low-order LB models can save computational costs without sacrificing simulation accuracy for 
gas flows with mixed Knudsen numbers. 

Discrete velocity sets can be obtained from several ways, see Refs. [2l)Il3"6"Il3"TIl4"T] . A direct method is utilizing 
the roots of Hermite polynomials [26] . In one-dimension, the discrete velocities £ a are just the roots of Hermite 
polynomials, and its corresponding weights are determined by: 

Wa = [»x"-(wr (15) 

For higher dimensions, the discrete velocity set can be constructed by using the "production" formulae 26 . 
Once the discrete velocity set is chosen, Eq.(|4]) can be discretized as 

^+€a-V/ a = --^(/ a -/^)+ 5a) (16) 

where f a — Wa J/l f ' , fa q = W "\,i£) an d 9a — "°Jt) ■ Therefore, the LB equation, i.e. Eg. ([To]) , is 
now obtained by discretizing Eq.(j4l) in the velocity space. 

2.2 Coupling scheme 

The key to success of a coupling scheme is appropriate bi-directional extraction and transfer of information 
at the interface or "hand-shaking "region. Since only LB models are used here, the extraction and transfer 
of information is in principle seamless. Lower order LB models, in their applicable capacity, can also retain 
non-equilibrium information which is required by the higher-order models. For instance, the D2Q16 model can 
already perform well for a range of Knudsen numbers [371 HH I2S1 131] • The D2Q9 model, which has been used 
for hydrodynamic simulations, may also capture some non-equilibrium effects [271 129 ) . This is very different 
from the continuum methods in kinetic-continuum hybrid models where non-equilibrium information is lost. 



1 We follow the conventional terminology for the LB models as first introduced in Ref. 1431 dubbed as DnQm model i.e. n 
dimensional model with m discrete velocities 



Interface 
High order LB model 



Low order LB model 

Figure 1: Schematic illustration of interface treatment where I is the grid on the interface while h and I represent 
the adjacent grids at the computational domains for high-order and low-order LB models respectively. 

To correctly transfer information across the interface between two LB models with different discrete velocities, 
the interface can be treated as a "virtual boundary" . Since a properly determined interface should be located 
at smooth regimes where lower-order models are valid, the relevant information can be obtained by using 
extrapolation and interpolation techniques. Firstly, the related macroscopic quantities can be calculated by 
interpolation, so that the equilibrium part of "boundary conditions" is able to be obtained by the Maxwell- 
Boltzmann distribution. For the non-equilibrium part of information, recall that not only higher-order but also 
lower-order LB models can produce accurate non-equilibrium information in the interface flow region. Moreover, 
the information provided by two models should be completely same on the "boundary" . Therefore, the non- 
equilibrium part of information on the "boundary" for the low-order and high-order LB models can be obtained 
via extrapolating information on the grids adjacent to the "boundary". It is interesting to note that similar 
techniques have been used to construct the no-slip boundary condition for continuum problems [44] [45] 06] . 

To illustrate the scheme clearly, it is convenient to discuss in details for a one-dimensional example. However, 
the same methodology can be generalized for multi-dimensional problems. The distribution function can be 
decomposed into its equilibrium (f^(I,t)) and non-equilibrium (/™ e? (/,£)) parts, i.e., 

f a (I,t) = f?{I,t) + fT(r,t), (17) 

where the letter I denotes an interface grid (see FigfT]). The velocity direction needs to be further classified 
since different discrete velocity sets are used across the interface. Hereinafter, the symbol + denotes the discrete 
velocities (see FigfT]) pointing to the lower-order LB model side, — to the higher-order LB model side. Firstly, the 
macroscopic quantities related to the equilibrium distribution can be obtained simply by the linear interpolation, 
i.e., 

Pi+Ph na . 

Pi = 2 — ' ( 8 ) 

ui = 2 ' ( 1Q ) 

where I and h are the interface neighboring grids (see Fig(TJ) . If necessary, the temperature can also be obtained 
similarly. With these quantities, the equilibrium distribution can be written as 



f a q + (I,t) « w a+Pl \ 1 + £ a+ • ui + - [(&,+ ■ ujY -uj + (T- !)(£+ -D)]\, (20) 



fj_(l, t) « w a -Pi \ 1 + £ a _ • Ul + - [(€„_ ■ Ui y - U j + (T- !)(£_ -D)}\. (21) 



Note, £ Q _, £ a + and ui Q _, w a + belong to two different discrete velocity sets. Based on the equilibrium 
distribution functions, the required information can be transferred cross the interface. Meanwhile, a first order 
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Figure 2: Schematic diagram of square lattices. 



extrapolation scheme is employed to supplement the information for the non-equilibrium part, i.e. 

f:i q (i,t)=f a+ (i,t)-r a i(i,t), 

f:Z q (I,t) = f a -(h,t)-f«1(h,t). 



(22) 
(23) 



Therefore, the general process of the present multiscale LB simulation starts from initialization to get all 
the necessary information, e.g., the velocity field by utilizing either the lower-order or higher-order model. The 
next step is to decompose the computational domain and determine the coupling interface by choosing an 
appropriate switching criterion. The final step is to implement the multiscale computation with lower-order 
models for continuum or near-continuum regime, and higher-order models for more rarefied regimes. Two 
models with different discrete velocity sets are coupled on the interface as described above. The second and 
third steps are repeated until the converged solutions are obtained. 

The determination of interface, i.e., choosing an appropriate switching criterion (also called 'breakdown 
parameters'), is important to any coupling/hybrid strategy. Several parameters have been proposed in literature, 
e.g., the local Knudsen number based on the local spatial gradients of hydrodynamic variables Kul = 4 -f- 
{4> is the interested flow quantity, typically density, temperature or pressure) [57], the 'B' parameters B = 
max{\Tij\ , \qi\} (jij is stress and %, heat flux[48j). However, different parameters may give significantly different 
values. So defining an appropriate switching criterion remains an interesting problem in itself [21 113). Here, we 
do not intend to investigate the switching criterion and will use what has been reported in the literature. In 
the next section, we will focus on numerical test of the present multiscale LB method. 

3 Numerical simulations and discussion 



3.1 Numerical scheme 

To solve Eq. (|16[) . various numerical schemes can be used. For instance, if the first-order upwind finite difference 
scheme is chosen, one can obtain the standard form of LB model, i.e., the stream-collision mechanism. However, 
for some high-order LB models, the discrete velocity points do not coincide with the lattice points. One has to 
choose a numerical scheme to break the tie between the time step and the lattice spacing used in the standard LB 



simulation [57]. As some discontinuities may occur at wall surface in the following simulations, we will employ 
the forward Euler method for time discretization and the second-order total variation diminishing (TVD) scheme 
for space discretization along the grid (see FigJ5]) for Eq. (fT51 l[28, 5 2] [531 [M]. According to the characteristics of 
problems, one can also choose any other appropriate numerical method to solve Eq. (|16[) . 

Let /"'! denote the distribution function value f a at the nth time step in the node (xi, yj) (see Figj2]), the 
scheme can be written as 
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■i a 



n+l,j 



J a,i 
£,ayOt 



£ ax St 



-p n ,3 _ pn,3 

" r a,i+\/1 • / a,i-l/2 



(24) 
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T n,] + l/2 _ pn,j-l/2 
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Ua.i - J a.i ) + 9aO t 



where 5 X and S y are the uniform grid spacing, and fit is the time step, £ ax and ^ ay denote the phase velocity 
component at the x and y coordinates. The outgoing and incoming fluxes in the node (i,j) (see Figj2j) are 

£,axfit 



•77" J 
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f'i-,3 
J a : i 
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where 
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e nj = 



Ja,i J a.i— 1 

f"J _ f n <3-l 
J a.i J a.i 



and the minmod flux limiter is 



*(0) = max [0, min(l, 0)] 



(25) 

(26) 
(27) 
(28) 

(29) 
(30) 

(31) 



3.2 Diffuse reflection boundary conditions 

Boundary treatment is of importance to correctly capture non-equilibrium effects, e.g. flow characteristics in 
the Knudsen layer. The simple diffuse reflection model, which was developed by Maxwell in 1879 [IH], has been 
proved to be sufficiently accurate for flows over a broad range of Knudsen numbers. The LB version of the 
Maxwellian model has also been developed |50j . Its specific numerical implementation on LB simulations has 
been discussed in Refs. |51ll52] . In this work, the Version 1 of boundary conditions in Ref. 52 will be employed. 
For convenience, we assume 



S « w a |l + £ a • u + \ [(ia ■ uf -u 2 + (T- 1)(£ - £>)] 



(32) 



i.e., f^ q = pS. As the discretization is conducted along a Cartesian coordinate system (see FigJ2]), the treatment 
of wall boundary can be described as 



fa.k — Pw,kS(T Wik ,u Wtk ) £ a ■ n > 



(33) 
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Figure 3: Schematic illustration of wall boundary treatment. 
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«a-n)<0 



"l/I* 



E l&» ■n|5(T w;fc ,«w,fc)' 

(i a -n)>0 



(34) 



where the subscript W denotes the computational nodes at the wall, pw,k denotes the density on the wall nodes 
k (see Fig|3]), Tw,ki the temperature, Uw,k, the velocity, n, the unity normal vector to the wall. Here, the 
distribution function in the ghost nodes are assumed to be identical to those on the corresponding wall nodes. 



3.3 Kramers' problem 

The classic Kramers' problem is often used to assess model capability in capturing the flow characteristics in 
the Knudsen layer (up to a few mean free paths away from the wall). In this problem, a gas fills the half-space 
(y > 0) bounded by a plate at y = 0. A constant shear rate is applied along the plate at y —t oo. With this 
special setup, one can investigate the nonlinear Knudsen layer in detail. To correctly predict this Knudsen 
layer, a kinetic method is required. However, for the flow region far from the plate wall, a continuum method 
is sufficient. Therefore, the problem is appropriate to test the coupling approach described in Sec. 12.21 

In the simulations, the plate is fixed at y — and a constant shear rate is applied at y — 200A(A denotes 
the mean free path). The Maxwellian diffuse reflection boundary condition is employed for the fixed wall. The 
D2Q36 LB model is used for the region near the plate (up to 10A from the wall) and the D2Q9 model for the 
other regiorQ The results in Fig0] show that the nonlinear velocity profile is captured well by the multiscale 
LB method. It indicates that the coupling process can effectively exchange information bi-directionally. 



3.4 Steady Couette flow 

With a simple geometrical configure, Couette flow represents many realistic shear dominant applications, e.g. 
reader heads of a hard-disc driver, micro turbines and gas bearings. Moreover, Couette flow is a theoretically 
well defined problem. Therefore, it is generally used as a benchmark problem. Particularly, its geometry is so 
simple that the coupling interface can be determined easily. As has been shown by FigJSl the 'B'profiles [48 
indicate that the flow regimes near the wall are highly rarefied and the discontinuities occur at the wall. So we 



2 For convenience, the multiscale LB model will be named according to the rule M-D2Qi-j where i and j denote the discrete 
velocities of lower-order and higher-order model respectively. 
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Figure 4: The velocity profile of Kramers' problem. The symbols are the data from Loyalka et al. [55] . where 
the BGK equation was directly solved. 







Kn=0.1 — i — Kn=0.3 - Kn=0.8 — — - 




Kn=0.2 « Kn=0.4 


0.9 










■ 


0.8 






0.7 




- 


0.6 


-'B 


.0- 


m0.5 


-*. 


30 ' 0D -O. Go n oBO° D '° ,.»- 


0.4 


^ 


*""*""*»*,«»l°°°°™'° aD °° a ° B °°°°°°°l..».-«--»-**'"'** X 


0.3 




^ X -"-^^ >t ^ M ^^ K _^^^^^^^^^'^'''' 


0.2 




***"- 


0.1 




- 



0.1 0.2 0.3 0.4 0,5 0.6 0.7 0.8 0.9 1 



Figure 5: The 'B'profile [35] for the Couette flows at different global Knudsen numbers. 



can use higher-order LB models in the near-wall regions and lower-order models in the middle. Therefore, the 
Couette flow can serve as a benchmark problem for testing the coupling strategy. 

In the following simulations, the lower-order LB model will be employed for 70% of the computational region 
in the middle and the higher-order model for the other regions adjacent to the walls. The upper and lower 
plates are set to be moving oppositely with the same velocity magnitude, and the diffuse boundary condition is 
used for gas/wall interactions. 

In Fig|6j it is clear that the D2Q9 model is unable to describe the Knudsen layer, which was also reported 
previously [27] [2Hl US]) while the M-D2Q9-36 model can obtain satisfactory results with the global Knudsen 
number up to 0.5. When the global Knudsen number is larger than 0.5, the multiscale method starts to deviate 
more from the linearized BGK (LBGK) results. This is not surprising since the Knudsen layers overlap and 
rarefaction effect becomes important for the whole flow domain (cf. FigJS]). Note, a typical Navier-Stokes and 
DSMC hybrid model usually become problematic when the Knudsen number is over 0.1, e.g., see Fig. 4 in Ref.[3]. 
To some extent, this indicates the advantage of coupling the kinetic-based LB models. 

As has been shown J5T1 [2"8I f2"Tl l2"9"] . various higher-order LB models can satisfy different requirement on model 
accuracy in terms of capturing high-order rarefaction effects. Therefore, it is possible to choose LB models with 
appropriate discrete velocity sets according to the requirements of model accuracy and computational cost. For 
instance, although the D2Q36 model is used for the regions near the wall in the above simulations, the D2Q16 
model may also be able to perform well for Knudsen numbers up to 0.4, see FigJT] Therefore, there is some 
flexibility in choosing various-order LB models for the present multiscale method. 





Figure 6: Nondimcnsional velocity profiles for planar Couette flows where the velocity is normalized by the 
velocity difference between the two plates. 





Figure 7: Nondimcnsional velocity profiles for the planar Couette flows where the velocity is normalized by the 
velocity difference between the two plates. 
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Figure 8: Nondimensional dynamics velocity profiles for Oscillatory Couette flows where the velocity is normal- 
ized by the velocity amplitude of oscillating plate. 
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Figure 9: Nondimensional dynamics velocity profiles for oscillatory Couette flows where the velocity is normal- 
ized by the velocity amplitude of oscillating plate. 

3.5 Oscillatory Couette flow 

The oscillatory Couette flows can mimic flows in many microfluidic devices containing oscillating parts. Its 
setup consists of a stationary plate at y — lo and a moving plate at y = which oscillates harmonically in the 
lateral direction with velocity u — u w sin(wi). This flow can be characterized by the Stokes number 



P = 



IpujL 2 



(35) 



which represents the balance between the unsteady and viscous effects. Similar to the steady case, 70% of the 
computational domain is computed with the lower-order LB model. The results will be compared to those of 
the variance- reduced (VR) particle simulations and the VR method is discussed in the Ref . |56j . 

FigJ5]shows that both D2Q9 and M-D2Q9-36 models are valid in the hydrodynamic regime when the Knudsen 
number is low (Kn=0.0178). When the Knudsen number increases and the flows are in the transition regime, 
Figs JS] and [TU] show that the D2Q9 model along becomes inappropriate while the M-D2Q9-36 model still performs 
well. This demonstrates that the present multiscale method can work well for the flows with various degree 
of rarefaction. The simulation results of the M-D2Q16-36 model as presented in FigJTT] further indicates the 
flexibility in choosing various-order LB models. 

It is also interesting to investigate the computational performance of the multiscale method. We will test the 
computational performance of the D2Q9, D2Q36 and multiscale models. For the M-D2Q9-36, 10% of the flow 
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Figure 10: Nondimensional dynamical velocity profiles for oscillatory Couette flows where the velocity is nor- 
malized by the velocity amplitude of the oscillating plate. 
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Figure 11: Nondimensional dynamics velocity profiles for oscillatory Couette flows where the velocity is nor- 
malized by the velocity amplitude of the oscillating plate. 

region is computed by the D2Q36 model and the rest is simulated by the D2Q9 model. The simulations are run 
on a four-core PC (Intel Core 2 QuadQ6600@2.4GHZ) without parallelization (i.e., only one core is utilized). 
The time required for each computational step is 0.114 ms for the M-D2Q9-36 model, 0.250 ms for D2Q36 model 
and 0.097 ms for D2Q9 model respectively. Therefore, the present multiscale approach can effectively reduce 
the computational costs for mixed-Kn flows. Similar to other hybrid methods, the performance of multiscale 
approach depends on how the computational domain is divided and calculated by the lower and higher order 
LB models. However, the LB framework can have some flexibility since various discrete velocity sets can be 
chosen to satisfy the requirement on model accuracy at the minimum computing cost. The details about how 
to choose appropriate LB models can be found in many references e.g., Refs. |26ll27l l28l l29l l3il l36l 137] 

4 Concluding remarks 

A multiscale LB method utilizing low-order and high-order LB models has been developed for gas flow simu- 
lation. As a hierarchy of LB models with various discrete velocity sets can be chosen, the multiscale method 
offers flexibility in designing coupling strategy to strike appropriate balance between model accuracy and com- 
putational efficiency. The present coupling process is simple by using interpolation and extrapolation processes. 
Therefore, the difficulty associated with kinetic-continuum hybrid models which couple two different methods 
becomes amenable. Furthermore, the present methodology can be extended to develop other kinetic-kinetic 
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hybrid models e.g. using discrete velocity models. 
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